Lyapunov Modes for a Nonequilibrium System with a Heat Flux 

Tooru Taniguchi 1 and Gary P. Morriss 2 

'-The Rockefeller University, 1230 York Avenue, New York, NY 10021, USA. 
2 School of Physics, University of New South Wales, Sydney, New South Wales 2052, Australia. 



Abstract 

We present the first numerical observation of Lyapunov modes (mode structure of Lyapunov vectors) in a 
system maintained in a nonequilibrium steady state. The modes show some similarities and some differences when 
compared with the results for equilibrium systems. The breaking of energy conservation removes a zero exponent 
and introduces a new mode. The transverse modes are only weakly altered but there are systematic changes to 
the longitudinal and momentum dependent modes. 



1. Introduction 

The difference of two trajectories starting from infinitesimally nearby initial conditions, which is called the 
Lyapunov vector, plays an essential role in the description of stability or instability in dynamical systems. The 
exponential rate of expansion or contraction of the absolute value of the Lyapunov vector is the Lyapunov exponent, 
and its positivity means that the system has a dynamical instability and is called chaotic. The Lyapunov vector 
is introduced in each independent direction of phase space, so in general we have to consider a set of Lyapunov 
exponents, called the Lyapunov spectrum, and the corresponding Lyapunov vectors in a high-dimensional chaotic 
system. An algorithm to calculate the Lyapunov exponents and vectors in many-body systems has been developed 
by Benettin et al. [1,2,3], and Shimada and Nagashima [4], and the behavior of Lyapunov vectors has been 
investigated from various points of view, for example, the conjugate pairing rule for Lyapunov spectra in some 
thermodynamic systems [5,6,7,8], and the localization behavior of Lyapunov vectors [9,10,11,12,13,14,15,16], etc. 

Recently, wavelike structure of Lyapunov vectors, called Lyapunov modes, and the corresponding stepwise 
structure of Lyapunov spectra has drawn attention in many-particle chaotic systems [17,18,19,20,21,22,23]. This 
structure appears in the region of Lyapunov exponents with small absolute values. Lyapunov exponents have the 
dimension of inverse time, so the structure in the small Lyapunov exponents is supposed to be a reflection of 
the slow and macroscopic behavior of many-particle systems. Since the first observation of Lyapunov modes in 
computer simulations of hard disk systems [18,24,25,26,27] they have been found in systems with soft potentials, 
independent of dimensionality and system geometry. It is now believed that the zero modes arise due to the 
presence of conserved quantities [28], and the spatial dependence of the non-zero modes is essentially analogous 
to higher order Fourier components. It has been shown that the breaking of conserved quantities, by for example 
a change from periodic to hard wall boundary conditions, leads to the absence of particular modes. It is also 
known that there are stationary modes and time-oscillating modes among the non-zero Lyapunov modes, and the 
time-oscillating period of the Lyapunov modes is connected to that of an equilibrium momentum auto-correlation 
function [26,28]. 

In this paper we consider the effect of breaking energy conservation by adding energy at one end of the system 
and removing it at the other. Although this leads to a flow of heat energy through the system, this is not intended 
to be a realistic model of low dimensional heat flow. Quasi-one-dimensional systems have been extensively used to 
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investigate Lyapunov modes and the localization or derealization of Lyapunov vectors for many-particle systems 
at equilibrium [25,26,15,16], and these systems can easily be adapted to maintain a heat flow. Using such a system, 
we investigate nonequilibrium effects of the Lyapunov modes and the momentum auto-correlation functions. 



2. Quasi-One-Dimensional Heat Model 

We consider a quasi-one-dimensional many-hard-disk system with a mechanism to insert energy at one end of 
the system and to remove energy from the other end. It consists of N hard disks of radius R and mass m in a two- 
dimensional rectangular region of length L x and width L y , with L y < 4R so that the disk order remains invariant. 
We choose periodic boundary conditions in the transverse y-direction, and in the longitudinal a;-direction we use 
a special boundary condition that allows the transfer of energy between the system and walls: 

p' x = -(l-e)p x -eVmk^Tsga(p x ) (1) 

Py = Py (2) 

where sgn(a:) = x/|a;|, Ub is Boltzmann's constant and e is a parameter (0 < e < 1). In a collision, (px,p y ) is 
the incoming momentum of the particle and (p' x ,p' y ) is the outgoing momentum. The parameter e controls the 
strength of the coupling to the "heat reservoir" of temperature T — T a (where a is either L or R). When e = 1 
the incoming momentum p x is completely replaced by the mean thermal momentum —y/rnksT sgn(p x ) of the 
reservoir and all information contained in the incoming momentum is lost. At e — the connection with the 
reservoir is removed and the collision process is hard wall. Here we consider values of e between and 1 so that 
the reservoir is coupled to the system but the loss of information is not complete. We choose the temperatures 
of the reservoirs independently, so the temperature of the left-hand side reservoir is Tl = 10, and the right-hand 
side reservoir is Tr — 1. A schematic illustration of this system is given in Fig. 1. For the numerical results that 
are shown in this paper, we use the values m = 1, N = 100, R = 1, L y = 272(1 + 10~ 6 ) and L x = 1.5NL y . 



Fig. 1. The quasi-onc-dimcnsional hard-disk system with a heat flux. The system width is narrow so that the disks remain in the 
same order (numbered 1, 2, ■ ■ ■ , N from left to right). When disks collide with walls at cither end of the system energy is transferred 
between the disks and walls depending upon the values of the temperatures Tl and Tr. 

Note that boundary condition (2) with periodic boundary conditions in the y-direction guarantees that the total 
momentum in the y-direction is conserved. The parameter e determines the strength of the coupling between the 
reservoir and the system, and specifies the fraction of incoming momentum p x that is retained after the collision. 
If Tl 7^ Tr and e / then an energy transfer (heat) occurs across the system, and a nonequilibrium steady state 
is established after a long time. 



3. Stepwise Structure of Lyapunov Spectra 

Figure 2 is a graph of the absolute value of the normalized Lyapunov spectrum {|A'"'|/A^} n with e = 0.5. Here 

the largest Lyapunov exponent is A^ 1 -* ~ 1.79. From this figure we see that the positive branch and the negative 

branch of the Lyapunov spectrum are not symmetric, namely A"' 7^ — a( 4JV_ -'+ 1 ) particularly for the exponents 

of smallest magnitude near N — 200. This is different from the equilibrium case (e = 0) where the dynamics 

is Hamiltonian with the symplectic structure guaranteeing the conjugate pairing rule for the Lyapunov spectra 
A G) = _ A (4iv-j+i) j 29L The system has 3 

zero-Lyapunov exponents, due to total momentum conservation and 
spatial translational invariance in the y— direction, and time-translational invariance. In the nonequilibrium case 
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Fig. 2. The Lyapunov spectrum for the quasi-onc-dimcnsional heat model consisting of 100 disks with e — 0.5. The reduced 
Lyapunov exponents ({ | A^ n ' | /A^ 1 -* } n ) with smallest magnitude are shown in the main plot. The inset is a plot of the positive 
exponent region of the full Lyapunov spectrum. The circles arc the smallest positive exponents, while the crosses are the absolute 
values of the conjugate negative exponents. The horizontal axis is the exponent number n, numbered from 1 (the largest) to 
4N = 400 (the smallest). 



Table 1 

Classification of the exponents in the Lyapunov spectrum for the quasi-one-dimcnsional system of 100 particles. The four zero 
exponents are 199. 200, 201 and 202 at equilibrium. 





Positive branch 


Negative branch 


Step number 


1-point step 


2-point step 


1-point step 


2-point step 


First 


198 


197, 196 


203 


204, 205 


Second 


195 


194, 193 


206 


207, 208 


Third 


192 


191, 190 


209 


210, 211 



(e 7^ 0) the total energy of the system is no longer conserved and this is the reason that A' 202 ' 7^ 0, although it is 
zero at equilibrium. 

In Refs. [25,26,28], it is shown that for the equivalent equilibrium system (that is, e = and hard- wall boundary 
conditions in the longitudinal direction and periodic boundary conditions in the transverse direction), the Lya- 
punov spectrum consists of 1-point steps and 2-point steps. The 1-point steps correspond to Lyapunov vectors 
containing transverse Lyapunov modes, and have their origin in total momentum conservation (and translational 
invariance) in the y— direction. On the other hand, the 2-point steps correspond to longitudinal and momentum 
proportional mode contributions. In Fig. 2 we can see that, away from equilibrium, the 1-point steps and 2-point 
steps remain (see Table 1 for meaning of 1-point and 2-point steps of the Lyapunov spectrum). The positive and 
negative branch 1-point steps are symmetric, that is A (n) « _a (4A, ~' i+1) , for n = 198, 195, 192, 189, 188, 185, • • ■. 
This may be because total momentum conservation and spatial translational invariance in the y— direction is 
preserved regardless of the value of e. On the other hand, the step-heights of the 2-point steps are clearly different 
in the positive and negative branches. Indeed, the asymmetry of the positive and negative exponents appears to 
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0.2 



Fig. 3. The local time average (5x j )t of the longitudinal spatial component of the Lyapunov vector corresponding to Lyapunov 
exponent \( 202 > as a function of the collision number rit and the normalized local time average of the particle position (xj)t/L x , 
for e = 0.5. 

be of the same magnitude for each 2-point step, and as well as for exponent 202 that has shifted from zero. 

There is no Gaussian thermostat [6,7] in this system so conjugate pairing of Lyapunov exponents would not 
be expected to hold. The magnitudes of the largest and smallest Lyapunov exponents appear to be almost equal, 
and there is only weak evidence for phase space dimensional contraction because dissipation as the origin of such 
a contraction can occur only at the two particles at the ends of the system. 



4. Mode Structure of Lyapunov Vectors 

For each Lyapunov exponent in the Lyapunov spectrum there is a corresponding Lyapunov vector which will in 
general have both coordinate and momentum components. For the exponents of smallest magnitude, the Lyapunov 
vectors often contain either coordinates of momentum components that vary sinusoidally with particle position 
forming a delocalized structure. We refer to such Lyapunov vectors containing delocalized structures as Lyapunov 
modes. 

4.1. Modification of zero-Lyapunov Modes 

The introduction of the nonequilibrium boundary condition e / leads to a change in the structure of the 
Lyapunov mode corresponding to exponent A" 02 '. The x— component of the new mode is shown in Fig. 3. The 
fact that A' 202 ) 7^ 0, is a purely nonequilibrium effect. Figure 3 shows the longitudinal component &xf 02 ^ of the 
Lyapunov vector <5r' 202 ^ as a function of the collision number rit and the normalized particle position for e = 0.5. 
Note that in this paper the Lyapunov vector 8T^ corresponding to the Lyapunov exponent A'™' is normalized 
so that |5r'™'| = 1. In order to obtain a clear graph, we take a local time average, indicated by the notation 
(• • -}t, over 8iV collisions using data just after collisions. Figure 3 shows an approximately steady structure with 
a slope which is steepest near the two ends of the system. This structure is quite different from the equilibrium 
case, where the longitudinal component Sxf N+2) is found numerically to be almost zero in both space and time. 



4 



8=0.0 x 8=0.5 8=0.9 * 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

< Xj>/L x 

Fig. 4. The transverse Lyapunov modes for the quasi-onc-dimcnsional heat model as shifted plots of the normalized transverse 
Lyapunov vector component {Sy^ ) / {8y^ ) as functions of {5Xj )/L x for e — (crosses), e — 0.5 (circles) and e — 0.9 (triangles). 
The first 1-point step (n — 198) is shifted so that zero corresponds to 4 on the vertical axis, the second 1-point step (n — 195) 
is shifted so that zero corresponds to 2 on the vertical axis and the third 1-point step (n — 192) is not shifted. Notice that the 
transverse modes arc largely independent of the parameter e. 

4.2. Transverse Lyapunov Modes 

Each 1-point step of the Lyapunov spectrum has an associated Lyapunov vector containing a transverse mode. 
The y— component 8y^ for the j-th disk of the n th Lyapunov vector has a sinusoidal dependence on the particle 
position Xj/L x at equilibrium [25,26]. We show that this structure is only slightly changed when the system is 
driven away from equilibrium, e 7^ 0. In Fig. 4 we present the time average of the transverse mode [Sy^) as 

functions of the average position (5Xj) / L x of the j-th disk for different values of e. The first (198), second 
(195) and third (192) transverse modes are clearly visible for the nonequilibrium cases. Note that the transverse 
Lyapunov modes in Fig. 4 show a weak e-dependence. In the nonequilibrium case e / 0, there is a deviation from 
a sinusoidal curve in the transverse mode, while in the equilibrium case e — 0, the mode structure can be fitted 
nicely by a sinusoidal curve. Fig. 4 shows that the amplitude of the transverse mode decreases slightly from the 
high temperature region (left side) to the low temperature region (right side). This change in amplitude in the 
transverse Lyapunov modes is a nonequilibrium effect . 

4.3. Longitudinal and Momentum- Proportional Lyapunov Modes 

Some of the structure of the longitudinal modes in the system remains as the system departs from equilibrium. 
There is a 7r/2 phase shift in time between the Sx components of modes in the same 2-point step, that is between 
{6Xj } and (5xj 196 '), see Fig. 5. As the momentum components are approximately the same functional form 
as the spatial components, there is the same phase shift between (Sp^f ) an d {Sp^ 1 ^}. The same behavior is 
observed for the components of the modes in the first 2-point step in the negative branch of the spectrum, {Sx j 204 ') 

and (5Xj } differ by a ir/2 phase shift in time, as do the contributions (<Sp^ 04 ^) and (Sp^® 5 ^). In the positive 
branch 2-point step the spatial contributions are larger than the momentum contributions, in the negative branch 
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Fig. 5. Contour plots of longitudinal Lyapunov modes for the quasi-one-dimensional heat model for e — 0.5. Exponents 196 and 197 
are the first 2-point step in the positive branch of the spectrum and 204 and 205 are the first 2-point step in the negative branch 
of the spectrum. 



the relative sizes of the contributions are reversed. 

The principle new feature is that nodal lines are curved, with those contributions to the positive branch of the 
spectrum having the center of the system reaching the node later than the two ends. This gives the appearance 
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Fig. 6. Contour plots of momentum proportional Lyapunov modes for the quasi-one-dimcnsional heat model for e — 0.5. Notice 
that in all of these figures the nodal lines that would be at (xj)t/L x — 0.5 for e — arc shifted to a smaller values. 



of a "forward" moving wave in the positive branch modes and a "backward" moving wave in the negative branch 
modes. It is less apparent in Fig. 5 that the time-oscillating period in the negative branch is larger than that for 
the positive branch. Also, the 5x components in the negative branch do not have nodes at the two ends but these 
are shifted within the boundaries. 

The momentum proportional components of the modes in Fig. 6 are also similar to those at equilibrium. 
The phase differences of n/2 and the straight nodal lines remain the same. However, the central nodal line at 
{xj)t/L x = 0.5 is now shifted towards the high temperature end of the system for both {Syj/p y j)t and (5p y j/p y j)t- 
As with the longitudinal modes, the time-oscillating period in the negative branch is larger than in the positive 
branch. The other modes that are not shown in Fig. 6 have a similar structure, for example (5xj/p X j) for 196 and 
197 have a similar structure to {Syj /p y j) but with more apparent noise. Also, (Sp x j/p x j) for 204 and 205 have a 
similar structure to (8p y j/p y j) but with more apparent noise. 



5. Time-Oscillating Periods of Lyapunov modes and Velocity Auto-Correlation 



The period T ac f (in numbers of collisions) of the auto-correlation function (acf) C x for the longitudinal com- 
ponent of the momentum increases, and the amplitude decreases as a function of e (see Table 2). The time 
oscillating period of the Lyapunov vectors for the fc-th longitudinal (and the momentum proportional) mode 
in the positive branch of the Lyapunov spectrum decreases slightly as a function of e, whereas in the negative 
branch the period increases slightly as a function of e. Such changes of T^~ k ^ are especially large for small fc's, 
which correspond to the ones close to the zero-Lyapunov exponents. The relation T^ 1 ' = 2T ac f is not correct away 
from equilibrium, while it is justified in equilibrium [26,28]. This does not contradict the proof of this relation as 
that required equilibrium and is based on the functional form of the Lyapunov modes, in particular that they are 
purely sinusoidal. 
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Fig. 7. The oscillating part of the normalized velocity auto-correlation function C x as a function of the collision number n t . Here 
e varies from to 0.9. The values of C x are the arithmetic averages of the auto-correlation for 11 disks in the center of the system, 
and is normalized so that CV I n — 1. 



Table 2 

Time-oscillating period of the momentum auto-correlation function and Lyapunov modes 



e 





0.5 


0.9 


Tacf 


8210 


8560 


8830 


T (3) 


5700 


5700 


5600 


T (2) 


8300 


8200 


8100 


T (l) 


16400 


16200 


14900 


T (-l) 


16400 


19600 


56400 


T (-2) 


8300 


9300 


12500 


T (-3) 


5700 


6300 


8100 



6. Conclusion and Remarks 



The changes in the structure of the Lyapunov modes for a system maintained away from equilibrium by a 
boundary imposed heat flow are described in detail. The breaking of conservation of energy is shown to introduce 
systematic effects, changing the number of zero exponents, and modifying the form of other modes, with the 
basic structure of transverse, longitudinal and momentum proportional modes remaining. This work is largely 
descriptive and the theoretical origin of the changes in modes is yet to be developed, as it is for the equilibrium 
case. However, it is clear that the modes are not only equilibrium properties of many-particle systems, but also 
fundamental for nonequilibrium steady states. 

The original purpose in introducing the nonequilibrium boundary conditions [i.e. Eqs. (1) and (2)] was to keep 
dynamical properties, like the momentum conservation in the transverse direction and deterministic property of 
orbit, etc., in a nonequilibrium dynamics while breaking energy conservation. The equilibrium case is recovered in 
the limit as e — > 0. Therefore this choice of the boundary conditions was not intended to provide a realistic model 
for the interaction between a heat reservoir and a particle system. For example, a consequence of these boundary 
conditions is that the momentum distribution function for a particle in contact with the reservoir is very different 
from a Gaussian. Indeed, we found that the incoming momentum is close to Gaussian but the outgoing momentum 
is strongly peaked around the mean momentum of the reservoir, with the width of its distribution determined by 
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the degree of mixing with the incoming momentum (that is, determined by the value of e). 
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